/* ----------------------------------------------------------------------------
 ex26.C
 mbwall 24mar96
 Copyright (c) 1995-1996  Massachusetts Institute of Technology

 The code for this example is adapted from an original implementation
 by Thomas Grueninger (from Uni Stuttgart, visiting scholar at MIT)

 DESCRIPTION:
 GAs are lousy at solving the travelling salesperson problem.  But here is an
 example of how to do it anyway.  In this case we use a steady-state genetic
 algorithm and give you the compile-time choice of partial match or edge
 recombination crossover.
 This one looks really good when you draw an entire population of individuals
 and watch them all evolve in real time.  It becomes very obvious what the
 genetic algorithm is doing and how well the speciation is working (or not).
 This implementation is by no means efficient, so please do not complain
 about all of the silly little inefficient aspects of the implementation.  But
 it does get the job done.
 ---------------------------------------------------------------------------- */
#include <math.h>
#include <ga/GASStateGA.h>
#include <ga/GAListGenome.h>
#include <ga/garandom.h>
#include <ga/std_stream.h>
#include <time.h>
#include "ants.h"
#include "utilities.h"
#include "InOut.h"
#include "TSP.h"
#include "timer.h"
#include "ls.h"
#include "mpi.h"
#include "solution.h"
#define cout STD_COUT
#define cerr STD_CERR
#define endl STD_ENDL
#define ostream STD_OSTREAM
#define ifstream STD_IFSTREAM

// Set this up for your favorite TSP.  The sample one is a contrived problem
// with the towns laid out in a grid (so it is easy to figure out what the
// shortest distance is, and there are many different paths with the same
// shortest path).  File format is that used by the TSPLIB problems.  You can
// grab more problems from
//
//
// Apologies for using fixed-length arrays.  But this is an example, not
// production code ;)

//#define TSP_FILE "tsp_rect_20.txt"


// You can use either edge recombination crossover or partial match crossover.
// Which one you select makes a HUGE difference in the performance of the
// genetic algorithm.  Only one of the two following lines should be commented.
//#define XOVER PMXover       // (Partial Match Crossover)
#define XOVER ERXover         // (Edge Recombination Crossover)
float Objective(GAGenome&);
int Mutator(GAGenome&, float);
void Initializer(GAGenome&);
float Comparator(const GAGenome&, const GAGenome&);
int ERXover(const GAGenome&, const GAGenome&, GAGenome*, GAGenome*);
int PMXover(const GAGenome&, const GAGenome&, GAGenome*, GAGenome*);
void ERXOneChild(const GAGenome&, const GAGenome&, GAGenome*);

static int * visit;
static char ** CM;
static char ** CM1;
static char ** CM2;
long int termination_condition_ga(void)
/*
 FUNCTION:       checks whether termination condition is met
 INPUT:          none
 OUTPUT:         0 if condition is not met, number neq 0 otherwise
 (SIDE)EFFECTS:  none
 */
{
	return ((n_tours >= max_tours) || (elapsed_time() >= max_time)
			|| (best_til_now == optimal));
}
GASteadyStateGA GA_Init(int pop_size, int generations, float mutation,
		float crossover)
/*
 * Iinit GA state
 */
{
	cout << "The Travelling Salesman Problem (TSP) program.\n" << endl;
	visit = new int[n];
	CM = new char*[n];
	CM1 = new char*[n];
	CM2 = new char*[n];
	for (int i = 0; i < n; i++) {
		CM[i] = new char[n];
		CM1[i] = new char[n];
		CM2[i] = new char[n];
	}
	// See if we've been given a seed to use (for testing purposes).  When you
	// specify a random seed, the evolution will be exactly the same each time
	// you use that seed number.
	//seed
	srand(time(NULL));
	unsigned int seed = rand();

	GAListGenome<int> genome(Objective);
	genome.initializer(::Initializer);
	genome.mutator(::Mutator);
	genome.comparator(::Comparator);
	genome.crossover(XOVER);

	GASteadyStateGA ga(genome);
	ga.minimize();
	ga.pReplacement(1.0);
	ga.populationSize(pop_size);
	ga.nGenerations(generations);
	ga.pMutation(mutation);
	ga.pCrossover(crossover);
	ga.selectScores(GAStatistics::AllScores);
	//ga.parameters(argc, argv);

	cout << "initializing...";
	cout.flush();
	ga.initialize(seed);
	cout << "evolving...\n";
	cout.flush();
	return ga;
}
void GA_end() {

}
Solution PGA(int myid, int pop_size, int generations, float mutation,
		float crossover)
/*
 * Main routine for Parallel GA
 * take input is id of process
 * output is the solution
 */
{
	// Each process init its own population ,GA state and variables

	GASteadyStateGA ga = GA_Init(n_ants, generations, mutation, crossover);//();

	Solution * new_list_solution;
	Solution * list_solution;
	list_solution = new Solution[pop_size + 1];
	new_list_solution = new Solution[pop_size + 1];
	// Each process run its own ga steps and after each 10 iteration exchange population
	best_til_now = 1000000000;
	while (!ga.done() && !termination_condition_ga()) {
		if (max_time < elapsed_time()) {
			break;
		}
		ga.step();
		if (ga.generation() % 10 == 0) {
			if (is_opened == 0) {
				cout << "in generation " << ga.generation()
						<< " the best in process rank : " << myid << " is :"
						<< ga.statistics().bestIndividual().score()
						<< "\t time :" << elapsed_time() << "\n";
				cout.flush();
			} else {
				out_stream << "in generation " << ga.generation()
						<< " the best in process rank : " << myid << " is :"
						<< ga.statistics().bestIndividual().score()
						<< "\t time :" << elapsed_time() << "\n";
				out_stream.flush();
			}
			GAPopulation pop = ga.population();
			// now update the solution list
			for (int i = 0; i < pop_size; i++) {
				GAGenome * gen = &pop.best(i);
				list_solution[i].cost = (*gen).score();
				GAListGenome<int> & genome = (GAListGenome<int> &) (*gen);
				if (genome.head()) {
					for (int j = 0; j < n; j++) {
						int xx = *genome.current();
						int yy = *genome.next();
						list_solution[i].tour[j] = xx;
					}
				}

			}
			// now let's all reduce the population for next generation
			// receive buffer

			MPI_Allreduce(list_solution, new_list_solution, pop_size,
					MPI_Solution, MPI_SolutionMin, MPI_COMM_WORLD);
			// now all process has the best population , put it back int to its real ga population
			for (int i = 0; i < pop_size; i++) {
				GAGenome* newgene = &pop.best(i);
				GAListGenome<int> &sis = (GAListGenome<int> &) *newgene;
				while (sis.head())
					sis.destroy();
				for (int j = 0; j < n; j++) {
					sis.insert(new_list_solution[i].tour[j]);
				}

			}

			ga.population(pop);
			best_til_now = ga.population().best(0).score();

		}
	}
	const GAGenome genome = ga.population().best(0);
	if (is_opened == 0) {
		cout << "the shortest path found is " << genome.score() << "\n";
		cout << "total time taken is : " << elapsed_time() << "\n";
		//cout << "this is the distance from the sequence\n";
		//cout << genome << "\n\n";
		//cout << ga.statistics() << "\n";
	} else {
		out_stream << "the shortest path found is " << genome.score() << "\n";
		out_stream << "total time taken is : " << elapsed_time() << "\n";
		//cout << "this is the distance from the sequence\n";
		//cout << genome << "\n\n";
		//cout << ga.statistics() << "\n";
	}
	GAGenome * gen = &ga.population().best(0);
	Solution solution;
	solution.cost = (*gen).score();

	GAListGenome<int> & genome_list = (GAListGenome<int> &) (*gen);
	if (genome_list.head()) {
		for (int j = 0; j < n; j++) {
			int xx = *genome_list.current();
			int yy = *genome_list.next();
			solution.tour[j] = xx;
		}

	}

	best_solution.cost = solution.cost;
	for (int i = 0; i < n; i++) {
		best_solution.tour[i] = solution.tour[i];
	}
	return solution;
}
int PGA_main(int argc, char* argv[]) {
	if (is_opened == 0) {
		cout << "The PGA has been chosen.\n" << endl;
	} else {
		out_stream << "The PGA has been chosen.\n" << endl;
	}
	int myid, numprocs;
	// Init MPI environment
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);

	construct_datatypes();
	PGA(myid, 50, 10000, 0.1, 1.0);
	MPI_Finalize();
	return 0;
}

int GA_main(int argc, char* argv[]) {
	if (is_opened == 0) {
		cout << "The GA has been chosen.\n" << endl;
	} else {
		out_stream << "The GA has been chosen.\n" << endl;
	}
	GASteadyStateGA ga = GA_Init(n_ants, 10000, 0.1, 1.0);
	best_til_now = 1000000000;
	while (!ga.done() && !termination_condition_ga()) {
		if (max_time < elapsed_time()) {
			break;
		}
		ga.step();
		if (ga.generation() % 10 == 0) {
			if (is_opened == 0) {
				cout << "in generation " << ga.generation() << " the best is :"
						<< ga.statistics().bestIndividual().score()
						<< "\t time :" << elapsed_time() << "\n";
				best_til_now = ga.population().best(0).score();
				cout.flush();
			} else {
				out_stream << "in generation " << ga.generation()
						<< " the best is :"
						<< ga.statistics().bestIndividual().score()
						<< "\t time :" << elapsed_time() << "\n";
				best_til_now = ga.population().best(0).score();
				out_stream.flush();
			}
		}
	}

	const GAGenome genome = ga.statistics().bestIndividual();
	if (is_opened == 0) {
		cout << "the shortest path found is " << genome.score() << "\n";
		cout << "total time taken is : " << elapsed_time() << "\n";
		//cout << "this is the distance from the sequence\n";
		//cout << genome << "\n\n";
		//cout << ga.statistics() << "\n";
	} else {
		out_stream << "the shortest path found is " << genome.score() << "\n";
		out_stream << "total time taken is : " << elapsed_time() << "\n";
		//out_stream << "this is the distance from the sequence\n";
		//out_stream << genome << "\n\n";
		//out_stream << ga.statistics() << "\n";
	}
	GAGenome * gen = &ga.population().best(0);

	best_solution.cost = (*gen).score();
	GAListGenome<int> & genome_list = (GAListGenome<int> &) (*gen);
	if (genome_list.head()) {
		for (int j = 0; j < n; j++) {
			int xx = *genome_list.current();
			int yy = *genome_list.next();
			best_solution.tour[j] = xx;
		}

	}

	return 0;
}

// Here are the genome operators that we want to use for this problem.
// Thanks to Jan Kees IJspeert for isolating an order-of-evaluation problem
// in the previous implementation of this function.
float Objective(GAGenome& g) {
	GAListGenome<int> & genome = (GAListGenome<int> &) g;
	float dist = 0;
	if (genome.head()) {
		for (int i = 0; i < n; i++) {
			int xx = *genome.current();
			int yy = *genome.next();
			dist += instance.distance[xx][yy];//DISTANCE[xx][yy];
		}
	}
	return dist;
}

void Initializer(GAGenome& g) {
	GAListGenome<int> &child = (GAListGenome<int> &) g;
	while (child.head())
		child.destroy(); // destroy any pre-existing list

	int i, town;

	for (i = 0; i < n; i++) {
		visit[i] = 0;
	}
	town = GARandomInt(0, n - 1);
	visit[town] = 1;
	child.insert(town, GAListBASE::HEAD); // the head node

	for (i = 1; i < n; i++) {
		do {
			town = GARandomInt(0, n - 1);
		} while (visit[town]);
		visit[town] = 1;
		child.insert(town);
	} // each subsequent node
	//free(visit);
}

int Mutator(GAGenome& g, float pmut) {
	GAListGenome<int> &child = (GAListGenome<int> &) g;
	register int n, i;
	if ((GARandomFloat() >= pmut) || (pmut <= 0))
		return 0;

	n = child.size();

	if (GARandomFloat() < 0.5) {
		child.swap(GARandomInt(0, n - 1), GARandomInt(0, n - 1)); // swap only one time
	} else {
		int nNodes = GARandomInt(1, ((int) (n / 2 - 1))); // displace nNodes
		child.warp(GARandomInt(0, n - 1)); // with or without
		GAList<int> TmpList; // inversion
		for (i = 0; i < nNodes; i++) {
			int *iptr = child.remove();
			TmpList.insert(*iptr, GAListBASE::AFTER);
			delete iptr;
			child.next();
		}
		int invert;
		child.warp(GARandomInt(0, n - nNodes));
		invert = (GARandomFloat() < 0.5) ? 0 : 1;
		if (invert)
			TmpList.head();
		else
			TmpList.tail();

		for (i = 0; i < nNodes; i++) {
			int *iptr = TmpList.remove();
			child.insert(*iptr, GAListBASE::AFTER);
			delete iptr;
			if (invert)
				TmpList.prev();
			else
				TmpList.next();
		}
	}
	child.head(); // set iterator to root node

	return (1);
}

int ERXover(const GAGenome& g1, const GAGenome& g2, GAGenome* c1, GAGenome* c2) {
	int nc = 0;
	if (c1) {
		ERXOneChild(g1, g2, c1);
		nc += 1;
	}
	if (c2) {
		ERXOneChild(g1, g2, c2);
		nc += 1;
	}
	return nc;
}

void ERXOneChild(const GAGenome& g1, const GAGenome& g2, GAGenome* c1) {
	GAListGenome<int> &mate1 = (GAListGenome<int> &) g1;
	GAListGenome<int> &mate2 = (GAListGenome<int> &) g2;
	GAListGenome<int> &sis = (GAListGenome<int> &) *c1;

	int i, j, k, t1, t2, town;

	for (i = 0; i < n; i++) {
		visit[i] = 0;
		for (j = 0; j < n; j++) {
			CM[i][j] = 0;
		}
	}
	while (sis.head())
		sis.destroy();

	// create connection matrix
	mate1.head();
	for (j = 0; j < n; j++) {
		t1 = *mate1.current();
		t2 = *mate1.next();
		CM[t1][t2] = 1;
		CM[t2][t1] = 1;
	}
	mate2.head();
	for (j = 0; j < n; j++) {
		t1 = *mate2.current();
		t2 = *mate2.next();
		CM[t1][t2] = 1;
		CM[t2][t1] = 1;
	}

	// select 1st town randomly
	town = GARandomInt(0, n - 1);
	visit[town] = 1;

	for (int i = 0; i < n; i++) {
		CM[town][i] = 0;
	}
	sis.insert(town); // the head node

	GAList<int> PossFollowList;
	GAList<int> FollowersList[5];
	while (PossFollowList.head())
		PossFollowList.destroy();
	for (k = 0; k < 5; k++) {
		while (FollowersList[k].head())
			FollowersList[k].destroy();
	}

	// select the following town with the minimal no of next folling towns
	int nPoss, nFollow;
	for (i = 1; i < n; i++) {
		nPoss = 0;
		for (j = 0; j < n; j++) { // no of poss. following towns
			if (CM[j][town]) {
				nPoss += 1;
				PossFollowList.insert(j);
			}
		}
		// nPoss = 0;
		if (nPoss == 0) {
			do {
				town = GARandomInt(0, n - 1);
			} while (visit[town]); // no follower
			visit[town] = 1;
			for (int i = 0; i < n; i++) {
				CM[town][i] = 0;
			}
			sis.insert(town);
		} else {
			PossFollowList.head();
			for (j = 0; j < nPoss; j++) {
				nFollow = 0;
				town = (*PossFollowList.current());
				for (k = 0; k < n; k++) {
					if (CM[k][town])
						nFollow++;
				}
				FollowersList[nFollow].insert(town);
				PossFollowList.next();
			}
			k = 0;
			while (FollowersList[k].size() == 0)
				k++;
			FollowersList[k].warp(GARandomInt(0, FollowersList[k].size()));
			town = (*FollowersList[k].current());
			visit[town] = 1;
			for (int i = 0; i < n; i++) {
				CM[town][i] = 0;
			}
			sis.insert(town);
		}
		while (PossFollowList.head())
			PossFollowList.destroy();
		for (k = 0; k < 5; k++) {
			while (FollowersList[k].head())
				FollowersList[k].destroy();
		}
	}
	sis.head(); // set iterator to head of list
	//free(CM);
	//free(visit);
}

int PMXover(const GAGenome& g1, const GAGenome& g2, GAGenome* c1, GAGenome* c2) {
	GAListGenome<int> &mom = (GAListGenome<int> &) g1;
	GAListGenome<int> &dad = (GAListGenome<int> &) g2;

	int a = GARandomInt(0, mom.size());
	int b = GARandomInt(0, dad.size());
	int h;
	if (b < a) {
		h = a;
		a = b;
		b = h;
	}

	int* index;
	int i, j, nc = 0;

	if (c1) {
		GAListGenome<int> &sis = (GAListGenome<int> &) *c1;
		sis.GAList<int>::copy(mom);
		GAListIter<int> diter(dad);
		index = diter.warp(a);
		for (i = a; i < b; i++, index = diter.next()) {
			if (*sis.head() == *index) {
				sis.swap(i, 0);
			} else {
				for (j = 1; (j < sis.size()) && (*sis.next() != *index); j++)
					;
				sis.swap(i, j); // no op if j>size
			}
		}
		sis.head(); // set iterator to head of list
		nc += 1;
	}
	if (c2) {
		GAListGenome<int> &sis = (GAListGenome<int> &) *c2;
		sis.GAList<int>::copy(mom);
		GAListIter<int> diter(dad);
		index = diter.warp(a);
		for (i = a; i < b; i++, index = diter.next()) {
			if (*sis.head() == *index) {
				sis.swap(i, 0);
			} else {
				for (j = 1; (j < sis.size()) && (*sis.next() != *index); j++)
					;
				sis.swap(i, j); // no op if j>size
			}
		}
		sis.head(); // set iterator to head of list
		nc += 1;
	}

	return nc;
}

float Comparator(const GAGenome& g1, const GAGenome& g2) {
	GAListGenome<int> &a = (GAListGenome<int> &) g1;
	GAListGenome<int> &b = (GAListGenome<int> &) g2;

	int i, j, t1, t2;
	float dist = n;

	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			CM1[i][j] = 0;
			CM2[i][j] = 0;
		}
	}

	// create connection matrix CM1
	a.head();
	for (i = 0; i < n; i++) {
		t1 = *a.current();
		t2 = *a.next();
		CM1[t1][t2] = 1;
		CM1[t2][t1] = 1;
	}
	// create connection matrix CM2
	b.head();
	for (i = 0; i < n; i++) {
		t1 = *b.current();
		t2 = *b.next();
		CM2[t1][t2] = 1;
		CM2[t2][t1] = 1;
	}
	//calc distance = how many edges are different
	for (i = 0; i < n; i++) {
		for (j = i; j < n; j++) {
			if (CM1[i][j] & CM2[i][j])
				dist--;
		}
	}
	free(CM1);
	free(CM2);
	return (dist);
}

//   Here we override the _write method for the List class.  This lets us see
// exactly what we want (the default _write method dumps out pointers to the
// data rather than the data contents).
//   This routine prints out the contents of each element of the list,
// separated by a space.  It does not put a newline at the end of the list.
//   Notice that you can override ANY function of a template class.  This is
// called "specialization" in C++ and it lets you tailor the behaviour of a
// template class to better fit the type.
template<> int GAListGenome<int>::write(ostream & os) const {
	int *cur, *head;
	GAListIter<int> tmpiter(*this);
	if ((head = tmpiter.head()) != 0) {
		os << *head << " ";
		for (cur = tmpiter.next(); cur && cur != head; cur = tmpiter.next())
			os << *cur << " ";
	}

	return os.fail() ? 1 : 0;
}
// force instantiations for compilers that do not do auto instantiation
// for some compilers (e.g. metrowerks) this must come after any
// specializations or you will get 'multiply-defined errors when you compile.
#if !defined(GALIB_USE_AUTO_INST)
#include <ga/GAList.C>
#include <ga/GAListGenome.C>
GALIB_INSTANTIATION_PREFIX GAList<int> ;
GALIB_INSTANTIATION_PREFIX GAListGenome<int> ;
#endif
